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Abstract. In meteorological and oceanological studies the classical approach for finding the numeri- 
cal solution of the regional model consists in formulating and solving the Cauchy-Dirichlet problem. 
The related boundary conditions are obtained by linear interpolation of data available on a coarse 
grid (global data), to the boundary of regional model. Errors, in boundary conditions, appearing 
owing to linear interpolation may lead to increasing errors in numerical solution during integration. 
The methods developed to reduce these errors deal with continuous dynamic assimilation of known 
global data available inside the regional domain. Essentially, this assimilation procedure performs 
a nudging of large-scale component of regional model solution to large-scale global data compo- 
nent by introducing the relaxation forcing terms into the regional model equations. As a result, 
the obtained solution is not a valid numerical solution of the original regional model. In this work 
we propose the optimization approach which is free from the above-mentioned shortcoming. The 
formulation of the joint problem of finding the regional model solution and data assimilation, as a 
PDE-constrained optimization problem, gives the possibility to obtain the exact numerical solution 
of the regional model. Three simple model examples (ODE Burgers equation, Rossby-Oboukhov 
equation, Korteweg-de Vries equation) were considered in this paper. The result of performed nu- 
merical experiments indicates that the optimization approach can significantly improve the precision 
of the sought numerical solution, even in the cases in which the solution of Cauchy-Dirichlet prob- 
lem is very sensitive to the errors in the boundary condition. 

1 Introduction 

Studying and modelling different physical processes frequently require to solve the Cauchy-Dirichlet 
problem. In geophysical investigations, for solving of the Cauchy-Dirichlet problem, numerical 
methods are usually applied. The discrete form of equations and initial and boundary values at the 
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points of the grid are traditionally used in these methods. The problem is generally solved by using 
a proper numerical scheme. However, both initial and boundary values, which are obtained from 

25 measurements or other model outputs, contain errors, which can reach 30% from the true values. 

For some problems, the solution can be very sensitive to these errors. At the same time some 
values of the sought solution at some points of the inner region, along with initial and boundary 
data, are often available. These additional data also have errors. The question arises: Is it possible to 
improve the accuracy of the solution of the Cauchy-Dirichlet problem by using these additional data 

30 (which may contain errors)? For some well-known equations of mathematical physics and dynamic 
meteorology we will show here that the use of additional information on the solution values within 
the integration region can noticeably improve the accuracy of the sought solution. 

An improvement of the solution accuracy of the Cauchy-Dirichlet problem is extremely important 
in meteorology because this is the problem of regional weather forecast. A distinction between two 

35 types of atmospheric forecast models needs to be done. These are global models, making forecast 
for the whole earth, and regional models, which produce a forecast for a limited region. One of basic 
differences between these model types is the grid resolution. Global models use a low resolution 
space grid and regional models operate on a more dense mesh. The reason for the existence of these 
two model types is, basically, computational. Even nowadays it is not possible in acceptable time 

40 period to integrate global models with the detailed physics and with the space resolution of regional 
models. It is worth noticing that global and regional models include different physical processes. 
Global models describe large-scale slow-time varying processes, with the time period of more than 
3 hours and the space scale larger than 60 km. Regional models can simulate the evolution of 
mesometeorological fast processes (small cyclones, storms, tornados), with the time period of less 

45 than 3 hours and the spatial extent between 2 and 60 km. 

These models are described mathematically by systems of nonlinear partial differential equations. 
To solve the system corresponding to a global model we need only the initial condition and bottom 
and upper boundary conditions, as the domain is the whole sphere. To get the solution of a system 
corresponding to a regional model we also need lateral boundary conditions. Generally, the lateral 

50 boundary conditions for a regional model are obtained from the global model. Global model data 
corresponding to the domain of a regional model lateral boundary are interpolated into a regional 
mesh and then are used for finding of the Cauchy-Dirichlet problem solution. 

We should notice that the lateral boundary conditions obtained from the global model do not 
provide any information on the structures of the scales smaller than the size of the mesh of the global 

55 model. In fact, global models do not distinguish a local meteorological phenomenon of characteristic 
length scale smaller than 30 km, because the space step is greater than this value. On the other 
hand, the parts of space and time spectrums of a regional model solution which correspond to large- 
scale structures and long-period processes are in worse accordance with observations than those of 
a global model. This is related to the fact that a regional model does not have the information about 
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60 the phenomena that occur outside its domain and therefore cannot describe with necessary precision 
(because the lateral boundary data are in low resolution discrete space and time grid) the impact of 
these phenomena on the evolution of processes inside the integration domain. 

For example, the weakly regular, low frequency oscillations of the sea surface temperature near 
the Peru coast, known as El Nino and La Nina events, have a large influence on the climate of 

65 Brazil. The circulation regimes over the Brazilian territory noticeably differ during these events. As 
a result, the precipitation rate in the South of Brazil during the La Nina periods is two times smaller 
than during the El Nino periods. Regional models driven, for example, by the Reanalysis data do not 
reproduce with confidence this distinction between the El Nino and La Nina regimes. The reason of 
this inconsistency is, probably, poor space and time data resolution on the regional model boundaries. 

70 Mathematically, it means that we are solving boundary value problems with significant errors on the 
boundaries and the solution is sensitive to these errors. 

The spectral nudging technique is one of methods proposed to use additional data from the inner 
domain in order to reduce boundary error influence and to improve a solution of initial-boundary 



value problem (Waldron et all. 1996 Storch et all. 2000 Kanamaru and Kanamitsu 2005). This 



75 method supposes incorporation of the largest internal modes of some meteorological variables from 
observations or from a driving model into the regional model solution. However the use of the 
spectral nudging technique requires inserting additional forcing terms into the evolution equations 
of a regional model. Hence the original model is corrupted and its new solution may not be close to 
the exact sought solution. 

80 In this work we propose a method which is free from the above mentioned shortcoming. Firstly, 
in Section 2 we give the formulation of the Cauchy-Dirichlet problem for a regional model and the 
general formulation of an optimization problem, and we show how the latter (with the accessibility 
of some additional conditions) can be applied for finding the solution of the former. In Section 3, for 
some simple equations, we show how the optimization approach can be used to obtain a solution of 

85 the initial boundary value problem. The equations considered here include the ordinary differential 
equation of Burgers, the one-dimensional, linearized Rossby-Oboukhov partial differential equation, 
and the partial differential equation of Korteweg - de Vries. We also compared the sensitivity of 
the solutions, obtained by the traditional and the new approach with respect to the errors in the 
boundary conditions. In Section 4 we discuss the obtained results. Appendix A provides more 

90 detailed description of the optimization procedure. 

2 Numerical solution of regional problems 

2. 1 The numerical forecast model 

We formulate here the Cauchy-Dirichlet problem as it is formulated for the regional weather forecast. 
To find a unique solution for the regional model equations we need the lateral boundary condition 
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for entire time interval on which the solution is seeking. We obtain these data from the solution of an 
outer model which, as a rule, is a global model. We don't know an exact solution of the outer model 
(if we did, our problem already would have been solved), but we can find its approximate solution 
in the discrete form: 

A t m = F d , 

^>{x,Q) =Y glohal {x), (1) 
%(x,t) \ xeB = %(x,t). 

Here A t is the evolution operator of the discrete model, W is the vector of the prognostic functions, 
Fd are discrete external forces, Y g i i, a i is the initial condition, B are upper, lower (and maybe lateral) 
95 boundaries of the global model, and tyf, is the boundary condition for the global model. Let 

will be the solution of (QJ. We suppose that this solution is sufficiently close to the solution of 
hypothetical ideal model which exactly describes the real atmospheric processes. 

The regional model is located in a closed area with boundary S and its discrete representation can 
be written in the same manner as for the global model 

S t {G} = F rd , 

G(x,0) = Y local (x), (2) 
G(x,t) \ xeS = G s (x,t), 

where St is the evolution operator of the discrete regional model, G is the vector of the prognostic 
functions and F r d are discrete external forces, Yi oca i is the initial condition, and G s is the boundary 

1 00 condition for the regional model. 

Now, we can solve numerically the initial boundary value problem Q using the solution ^ s i of the 
global model ([T| for getting the boundary conditions. The traditional approach to seek the solution of 
Q consists in the interpolation of required data from ^ s i on the regional grid for forming the initial 
and boundary conditions and then applying any numerical method to solve the model equations. But, 

1 05 as we have mentioned above, the boundary conditions contain errors which can strongly corrupt the 
solution. 

The use of the information obtained from the solution of the global model equations fy s i on the 
values inside the area of the regional model integration for all available time moments can help to 
overcome this difficulty. 

1 1 2.2 Formulation of the optimization problem 

We assume that the solution of regional model Q cannot strongly differ from global model solution 
^si. Our objective is to find G s i on the fine (regional) mesh that satisfies the initial values and the 
equations of the regional model The solution G s i also have to be as close as possible to the 
ty s i inside the regional model domain. We can formulate this aim as an optimization problem in the 
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following manner 

Minimize d(G, ^f st ) 
subject to (s.t.) S t {G} - F rd = 0, (3) 
G(x, 0) - Yi oca i(x) = 0. 

Here d(:, :) is the objective function, which represents the distance between G and \I> s j vectors. In 
other words, we want to minimize the distance between the regional and the global model solutions 
under conditions that the regional model equations and the initial condition are satisfied. Note here 
that, when we talk about the solution G, we bear in mind the vector G = G(x, t) in a grid space of 

1 1 5 (x,t) - coordinates, where t corresponds to discrete time points of the regional model integration and 
x to discrete mesh points in the regional model area. 

There are different approaches to solve the optimization problem ([3j. For example, we can apply 
it to each time level, as it is done in finite-differences methods: For the given initial conditions 
(to = 0), the solution is obtained at t\. Then, considering the solution at t\ as the initial condition, 

1 20 the solution is found at t% and so on. On the other hand, we can consider the domain of definition 
of the problem on the mesh including all time levels. In this case all available information from the 
global model for the period of integration will be taken into consideration simultaneously. We use 
here this approach. 



3 Examples of application of the optimization method to the problems of small dimension 

125 3.1 Burgers' equation 

We shall demonstrate an application of the optimization method to problems of small dimension. 
As a first example we consider the supersensitive boundary value problem for the Burgers' equation 
dBohel[T996| : 

ex" — —xx', 

(4) 

x(0) = -l,x(T) = 1. 

To get the analytical solution of this equation it is enough to integrate the left and right sides of the 
equation Q over t: 

ex' 



x 2 + c 



2 

Then, after rewriting the foregoing formula as 

dx dt 
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and integrating the left side over x and the right one over t, we can write the Burgers' equation 
solution as 



-c I 1 + exp 



1 — exp 



for c < 0, 



/5tarx ^ t + £ 2 ^ V foroO. 

The graph of the analytical solution for equation Q with the boundary condition x(0) — — 1, x(l) = 
1 and e = 0.05 is presented in Figure[T| Numerical solution of this problem applying, for example, 
the shooting method coincides with great accuracy with the analytical one. 
[Insert figure [T] here.] 

130 Let us choose several points from the analytical solution and slightly perturb them (till 5% from 
its real value). This procedure models the boundary and inner domain data containing errors. Using 
these perturbed data for the boundary points x(0), x(l) we solve the boundary problem Q applying 
the shooting method. For the step At = 0.01 the numerical solution is presented in Figure|2]by red 
line. 

135 [Insert figure[2]here.] 

As one can see, small perturbations of the boundary condition result in large errors in the solution. 
Now, let us solve the optimization problem ([3]) using additional information about the perturbed 
solution on the points inside the domain. Using the same discretization as when solving |4} by the 
traditional shooting method, with the same At = 0.01 we find the solutions which are exhibited in 

140 Figure3. 

Figure [3] shows the solutions of the optimization problem for the three cases when we used 3, 4 
and 5 points of the perturbed solution, respectively. One can see that all these solutions correspond 
better to the analytical solution, than the solution obtained by the traditional approach. In the case of 
using the three additional inner points, the optimization solution and the analytical solution nearly 
145 coincide. This example shows that there are dynamical systems in which small perturbations (< 
5%) in the boundary conditions can lead to great errors in the final solution. However, if some 
additional information exists, it can be used to improve significantly the sought solution, applying to 
the considered problem the optimization theory methods. 
[Insert figure[3]here.] 

150 3.2 Rossby-Oboukhov one-dimensional equation 

Here we consider the one-dimensional linear partial differential equation obtained during the lin- 
earization procedure of the Rossby-Oboukhov equation ( Obouk hov] 1 1949 1: 



— (— --^U + /3— + U— -0 (5) 
dt \ dx 2 Iq J dx dx 3 
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Here if) is the stream function, f = 10~ 4 s _1 is mean value of Coriolis parameter, (3 = df /dy = 
1.6 • 10~ n s _1 m _1 is mean value of meridional gradient of Coriolis parameter, Iq — C0//0 = 3 ■ 10 6 
m is the Oboukhov scale, cq is the sound velocity, and U is the zonal wind, which is variated between 
and 30 m/s. 

For the periodic boundary condition 

iP{0,t)=i/;(L,t), (6) 

where L is the size of integration area (we shall use L — 3 ■ 10 7 m), the solution of equation ^ can 
be written as 

N 

t) = ^ A n sin[fc n (x - c n t) + 4> n ], (7) 

n=l 

where 

C-V-gM. w 
kl + l/l 2 

™ = ^zT' 

1 55 and A n and <f> n are defined by the initial condition (Rossby 1939 ). 

For finding the numerical solution it is convenient to rewrite the equation in nondimensional form. 
We choose the following scales 5 = 6- 10 6 m, T = S/V — 6 ■ 10 5 s, V = 10 m/s. The dependent 
and independent nondimensional variables are defined as follows: 

t=t i>=^- (9) 

, and equation |5]l may be written in nondimensional form as 

d ( d 2 _ ^-\7 + g9ip + u d 3 ij_ Q 
dt \ dx 2 b 2 J dx dx 3 

where § = £ = 2, /3 =^ = 57.6 and f/ =^e [0, 3]. 

For finite-difference discretization we use unconditionally stable scheme with truncation error 
0(Ax 2 , At 2 ) given by 

1 A^+ 1 1 -2^+ 1 +^+ 1 $n-2#+$L x 



At \ Ax 2 Ax 2 b 2 

2 I 2Ax 2Aa; y 



2Aa; 3 2Ax 3 



= 



At first, we generate a specific analytical solution (|7| containing 85 modes. Figure |4] shows this 
solution for t = h. 
[Insert figure [4] here.] 
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1 60 Note that in Figures we use dimensional values, but in the numerical computations all variables 
are nondimensional. 



As a local model we consider the equation (10 1, with initial and boundary conditions, defined 
over [a, b] £ [0, L) (closed interval smaller than entire domain). For initial and boundary conditions 
we will use our specific form (see Fig. [4]) of the global model solution ([7J. For simulating really 

1 65 encountered problems of atmospheric modelling we took the analytical solution in the points of the 
the coarse grid with space step Air = 200 km and time step At = 2 hour, and perturb its values 
randomly in such a way that a perturbation can reach till 30% of its exact value. 

Primarily we find the solution of the forward Cauchy-Dirichlet problem for our local model. As 
local domain we take the interval [1.8 • 10 7 , 2.4- 10 7 ]m (6000 kilometers) inside the global model in- 

1 70 terval [0,3- 10 7 ]m. For the first reference experiment the initial and boundary conditions were taken 
from the exact analytical solution Q. Requiring that the numerical solution with the exact initial 
and boundary condition have to nearly coincide with the analytical solution, we took as the maximal 
possible values of space and the time steps Ax = 10 kilometers and At = 200 seconds respec- 
tively. To obtain the solution of the forward problem described above for 96 hours it is demanded 

1 75 0.8 second of CPU time. 
[Insert figure|5]here.] 

In the next experiments the perturbed analytical solution on the coarse grid was used for formation 
of the boundary conditions. Figure [5] shows the exact solution and the solution with perturbed 
boundary conditions of the forward problem. 
180 [Insert figure [6] here.] 

One can see that for the forward problem the numerical solution rather quickly diverges from the 
analytical solution due to the errors of the linear interpolation procedure in the border points. After 
96 hours perturbed boundary solution has very faint resemblance with the analytical one. 

In the second series of experiments we applied the optimization method to the local model. For 
1 85 this calculation we used all available global data (perturbed analytical solution on coarse grid) in the 
inner local domain. Figure [6] shows the results of calculations. It is important to note, that, as in 
previous experiments, we chose the space Ax =100 km and the time At = 3600 sec steps to be as 
large as possible. 

[Insert figure [7] here.] 

190 In the optimization approach one can use larger steps than in the forward problem, because we 
do not deal with a time-evolution problem, and there is no accumulation of numerical errors at 
each time step. In Figure [7] one can see that for space and time steps many times larger than 
those that were used in the forward problem (Ax = 100 kilometers and At = 3600 seconds) we 
have far better agreement with the analytical solution. The decrease of time and/or space step in 

1 95 the optimization approach does not appreciably improve the solution, only the smallest details are 
slighly better reproduced. 
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As the local mesh in the optimization approach is rather coarse, the average CPU time needing 
for solving the problem is smoller than for forward problem (0.5 seconds) in spite of greater com- 
putational complexity. Also we have to pay attention on a very weak sensitivity to the global data 
200 perturbations. The impact of every individual perturbation on the optimization solution is very small. 
For example, increasing errors in randomly perturbed global data up to 60% does not strongly affect 
the solution as one can see in Figure [8] 

[Insert figure [8] here.] 

These experiments clearly show that the use of additional global data inside the local domain, 
205 even perturbed, can significantly improve the solution of the model. The optimization approach, 
albeit seeking the solution on entire space-time grid, is not very numerically expensive, because the 
space, and, principally, the time steps can be chosen much larger than in the forward problem. 

3.3 Korteweg-de Vries equation 

Now we consider the Korteweg-de Vries equation: 

u t + 6uu x + u xxx = 0, (11) 

and find its numerical solution applying both the traditional forward method and the optimization 
210 method. 

The equation ( fTT| ) has the exact solution ( jKorteweg and de Vries|[T895] l,( |Grimshaw|[2004[ ) 

u = 6 + acn 2 (7(x- Vt)\m), (12) 

where cn(x|m) is the Jacobi elliptic function, m (0 < m < 1) is the module of elliptic function, 
a = 2mj 2 and V = 6&+4(2m— l)j 2 . For the case when m — > 1 we will have cn(x\m) — > sech(x) 
and the solution ( p~2] > will have the form 

u = b + ascch 2 (j(x~ Vt)), (13) 

with V = 66 + 2a,and a = 2j 2 which describes a one-dimensional soliton . 



For the finite-difference discretization we use the following implicit numerical scheme (Furihata 



1999[ ), which possesses the properties of total energy and mass conservation: 



j + l 3 



- u 



At 



iV-i + («? +1 ) 2 -(«ti) a ) + 2^ 




I J+l i 3 \ i / J+l i 3 \ U i-2 ' U i-2 

-W+i + u i+i) + K-i + K-i) 

As a numerical scheme is non-linear for obtaining of numerical solution of a forward problem we 
use the Newton method. 
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Following the steps of previos example, we choose as a reference model the analitical solution 
215 with cofficient b = a, 7 = 2 and module of elliptic function m = 0.995. The Figure [9] represents 
the solution on domain [—5, 15] at time t = 0. 
[Insert figure [9] here.] 

As a local model we consider the equation (Hi defined over closed interval [0, 10]. To get a 
good accordance between a numerical solution of forward problem (with exact initial and boundary 
220 conditions) and an anlytical solution on the time interval < t < 1 we took as the maximum 
possible values of space and the time steps Ax = 0.02 and At = 0.0002 respectively. The CPU 
time required to find the solution is 15 seconds. 



[Insert figure 10 here.] 

As a global model solution we take the analytical solution, intoduced above as the reference 
225 model, with space step Ax — 0.5 and time step At = 0.005 and perturb its values till 10% from the 
exact ones. Using these perturbed global model solution to form boundary condition for local model 
we will solve initially the forward problem. 



[Insert figure 11 here.] 

In Figure [TT] one can see that the forward problem with small errors in the boundary conditions 
230 can produce inacceptable numerical solution. 

Now, we apply the optimization method to the local model using additional available for inner 



points information from global model solution. Figure 12 shows the solution of the optimization 
problem. As in previous experiment we chose the space (Ax = 0.2) and the time (At = 0.002) 
steps as large as possible. 
235 [Insert figure [T2|here.] 

It can be clearly seen the advantage of the optimization approach in this case. The difference 
between numerical and exact solution is very small in comparison with the forward problem calcu- 
lations. 



4 Discussions and conclusions 

240 The results of the numerical experiments presented here indicate that the optimization approach can 
significantly improve the precision of the sought numerical solution of the regional model when the 
boundary values have errors but the information on the behaviour of the sought solution in a number 
of inner points is available. Even in the cases in which the solution of Cauchy-Dirichlet problem 
is very sensitive to the errors in the boundary condition the use of optimization approach gives the 

245 possibility to construct the solution which are close to the analytical solution. We also made exper- 
iments with two-dimentional nonlinear Rossby-Oboukhov equations. The preliminary results also 
demonstrate that the use of the optimization approach significantly improve the numerical solution 
of boundary problem with errors in the boundary values. At the present time we are developing 
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economic algorithm that can be applied to any great computational problem. 



250 Appendix A 

Formulation of optimization problem 

We can express our problem as the following non-linear optimization problem with equality con- 
straints 



Minimize h\\u — V\\ 2 



(Al) 



s.t. h(u) = 0, 

where V represents the global data on the regional mesh and h(u) = \h\{u) h 2 (u) . . . h m (u)) T is 
a vector of the discretizations of regional equations at each point of space-time mesh of the regional 
255 model. 

A usual optimization technique to solve this kind of problem is to apply Newton's iteration method 
to the system of nonlinear equations arising from first-order necessary conditions, known as Karush- 



the problem (Al 



Kuhn-Tucker (KKT) conditions (Nocedal and Wright 1999). For instance, the KKT conditions for 
i are 

(ti - V) + h'{u) T \ = 
h(u) = 0, 



(A2) 



where 



h'(u) = 



( Vh^uf \ 
Vh 2 (u) T 



\Vh m (u) T / 

is the Jacobian matrix of the constraint function and A represent the vector of Lagrange coefficients. 



Each step of the Newton iteration associated with system ( |A2| > is defined as follows: 



J(u k , A fe ) 



Uk - V + h'(u k ) T \ k 
h(u k ) 




(A3) 



where J represents the Jacobian matrix of the system ( A2 1 



J(u,A) 



V 



h'{u) 



(A4) 



/ 



and W 2 hi(u), i = 1, . . . , m, are Hessian matrices of the constraints. 
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The Jacobian matrix ( A4 1 is a saddle point matrix and there are many methods that can be applied 
with assosiated linear system (Benzi, Golub and Liesen 2005 | l. But first, note that the derivative 
of the discretization vector of the constraints is described by a sparse matrix with block-diagonal 
260 structure, because this derivative is calculated at every point of time-space mesh. On the other hand, 



the first part K^ 2 hi (u) of Jacobian matrix ( A4 1 of the system ( A2 1 includes the calculation of 



the Hessians of the constraints discretization vector, which is a computationally expensive procedure 
because calculations have to be made for every Newton iteration. The resulting matrix is generally 
dense. To accelerate the calculations and reduce the consuming of the computer memory we use the 
265 following Jacobian approximation: 



B k = 



(A5) 



/ h'(u k ) T 
h'(u k ) 

The use of the approximation Bk for J(iik, A&) strongly simplifies the procedure of finding the 
solution of the linear system ( |A3| > ( [Benzi, Golub and Liesen||2005| l. 
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Figure Captions 



Fig. 1. Analytical solution of the Burgers' equation for e = 0.05 and the boundary condition: 

x(0) = -l,a;(l) = 1. 

Fig. 2. Forward problem solution with 5% perturbed boundary condition of Burgers' equation in comparison 
with exact analytical solution. 

Fig. 3. Optimization problem solution of Burgers' equation for three cases: 3 reference points; 4 reference 
points; and 5 reference points. 

Fig. 4. Analytical solution of Rossby-Oboukhov equation on periodic domain [0, L] at t = h with 85 modes. 

Fig. 5. Forward problem solution with exact boundary condition of Rossby-Oboukhov equation at 96 hours. 
Ax = 10 km, At = 200 sec. CPU time = 0.8 sec. 

Fig. 6. Forward problem solution with 30% perturbed boundary condition of Rossby-Oboukhov equation at 
t = 24, 48, 96 hours. Ax = 10 km, At = 200 sec. CPU time = 0.8 sec. 

Fig. 7. Optimization problem solution of Rossby-Oboukhov equation with 30% perturbed global solution at 
t = 24, 48, 96 hours. Ax = 100 km, At = 3600 sec. Average CPU time = 0.5 sec. 

Fig. 8. Optimization problem solution of Rossby-Oboukhov equation with 60% perturbed global solution at 
96 hours. Ax = 100 km, At = 3600 sec. Average CPU time = 0.5 sec. 



Fig. 9. Analytical solution of KdV equation as the cnoidal wave 1 12 i on the domain [—5, 15] at t = with 
7 = 2, b = a and module m — 0.995. 

Fig. 10. Forward problem solution with exact boundary condition of KdV equation at t — 1. Ax = 0.02, 
At = 0.0002. CPU time = 15 sec. 

Fig. 11. Forward problem solution with 10% perturbed boundary condition of KdV equation at time t = 0.2, 
0.5, 1. Ax = 0.02, At = 0.0002. CPU time = 23 sec. 

Fig. 12. Optimization problem solution of KdV equation with 10% perturbed global solution at time t — 0.2, 
0.5, 1. Ax = 0.2, At = 0.002. Average CPU time = 21 sec. 
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Fig. 1. Analytical solution of the Burgers' equation for e = 0.05 and the boundary condition: 

a;(0) = -l,x(l) = 1. 




Fig. 2. Forward problem solution with 5% perturbed boundary condition of Burgers' equation in comparison 
with exact analytical solution. 
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Fig. 3. Optimization problem solution of Burgers' equation for three cases: 3 reference points; 4 reference 
points; and 5 reference points. 
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Fig. 4. Analytical solution of Rossby-Oboukhov equation on periodic domain [0, L] at t = h with 85 modes. 
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Fig. 5. Forward problem solution with exact boundary condition of Rossby-Oboukhov equation at 96 hours. 
Ax = 10 km, At = 200 sec. CPU time = 0.8 sec. 
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6. Forward problem solution with 30% perturbed boundary condition of Rossby-Oboukhov equation at 
24, 48, 96 hours. Ax = 10 km, At = 200 sec. CPU time = 0.8 sec. 
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Fig. 7. Optimization problem solution of Rossby-Oboukhov equation with 30% perturbed global solution at 
t = 24, 48, 96 hours. Ax = 100 km, At = 3600 sec. Average CPU time = 0.5 sec. 
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Fig. 8. Optimization problem solution of Rossby-Oboukhov equation with 60% perturbed global solution at 
96 hours. Ax = 100 km, At = 3600 sec. Average CPU time = 0.5 sec. 




Fig. 9. Analytical solution of KdV equation as the cnoidal wave 1 12 i on the domain [—5, 15] at t — with 
7 = 2, b = a and module m — 0.995. 
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Fig. 10. Forward problem solution with exact boundary condition of KdV equation at t — 1. Aa; = 0.02, 
At = 0.0002. CPU time = 15 sec. 
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Fig. 11. Forward problem solution with 10% perturbed boundary condition of KdV equation at time t 
0.5, 1. Ax = 0.02, At = 0.0002. CPU time = 23 sec. 
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Fig. 12. Optimization problem solution of KdV equation with 10% perturbed global solution at time t 
0.5, 1. Ax = 0.2, At = 0.002. Average CPU time = 21 sec. 
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